Preparing the data

Loading the packages and data directories

options(width=108)
out=F
# The list of subjects, the order of conditions, and the thresholds are derived from Subjects.xlsx
library(xlsx)
library(ggplot2)
library(plyr)
library(reshape)
library(matrixStats)
library(gridExtra)
library(lme4)
library(lmerTest)
library(BayesFactor)
#library(splines)
db <- '/home/egor/Dropbox/' # on Linux
db <- '/Users/Egor/Dropbox/' # Windows
# db <- '~/Dropbox/' # on Mac
source(paste(db, 'Prog/R/myFunctions/pvalfn.R', sep=''))

Plot variables

# theme for plotting:
alpha <- .7
themefy <- function(p) {
    p <- p + theme_bw() + 
         theme(panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(),
            axis.text=element_text(size=8), axis.title=element_text(size=9),
            legend.text=element_text(size=8), legend.title=element_text(size=9),
            legend.key = element_blank(), legend.margin=margin(t=-.04, unit='in'),
            legend.background = element_rect(fill='transparent'),
            plot.title=element_text(face='bold'))
}
cc <- c('#e31a1c','#fdbf6f','#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99')
xLab <- 'Target Peak Time (s)'
yLab <- 'Log Contrast Threshold'
yLab2 <- 'Normalized Contrast Threshold'
colLab <- 'Mask Velocity\nBandwidth'
colLab2 <- 'Mask Velocity Bandwidth'
colLabType <- 'Mask Type'
yLim <- c(-1.5, -0.8)
yLimNorm <- c(-1.3, -0.7)
dodge <- position_dodge(width=0.1)
sumFn <- function(ss, subjStr='subj', xStr='targTpeak', grpStr='maskV'){
    sumSubj <- ddply(ss, c(subjStr, xStr, grpStr), summarise,
                     mnS=mean(threshMean), se=sd(threshMean)/sqrt(length(threshMean))) #, .drop=F)
    sumSubjMn <- ddply(ss, c(subjStr), summarise, mnStot=mean(threshMean)) # total mean across conditions per subj
    sumSubj <- merge(sumSubj, sumSubjMn)
    sumSubj$normS <- - sumSubj$mnS / sumSubj$mnStot # normalized subject mean
    sumSubj$seNorm <- NA
    # sumSubj$mnS[is.na(sumSubj$mnS)] <- 0
    sumGrp <- ddply(sumSubj, c(xStr, grpStr), summarise,
                  mn=mean(mnS), se=sd(mnS)/sqrt(length(mnS)),
                  norm=mean(normS), seNorm=sd(normS)/sqrt(length(normS)))
    sumGrp$subj <- 'average'
    sumSubj <- rename(sumSubj, c(mnS='mn',normS='norm'))
    sumComb <- rbind(sumGrp, subset(sumSubj, select=-mnStot))
    sumComb$se[is.na(sumComb$se)] <- 0
    sumComb
}
plotAve <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='maskV', 
                    xlab=xLab, ylab=yLab, collab=colLab, yStr='mn', seStr='se'){
    pss$yMin <- pss[,yStr] - pss[,seStr]
    pss$yMax <- pss[,yStr] + pss[,seStr]
    p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
                            ymin='yMin', ymax='yMax')) +
        geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
        scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
        geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
        labs(x=xlab, y=ylab, colour=collab) + #ylim(0,1) + 
        guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
    p <- themefy(p)
}
plotIndiv <- function(pss, subjStr='subj', xStr='targTpeak', grpStr='maskV', 
                    xlab=xLab, ylab=yLab, collab=colLab2, yStr='mn', seStr='se'){
    pss$yMin <- pss[,yStr] - pss[,seStr]
    pss$yMax <- pss[,yStr] + pss[,seStr]
    p <- ggplot(pss, aes_string(x=xStr, y=yStr, colour=grpStr, group=grpStr,
                            ymin='yMin', ymax='yMax')) + 
        facet_wrap( ~ subj, ncol=4) +
        geom_point(position=dodge, size=1, alpha=alpha) + geom_line(position=dodge, alpha=alpha) +
        geom_linerange(position=dodge, show.legend=F, alpha=alpha) +
        scale_x_continuous(breaks=c(.5,1,1.5), labels=c('0.5','1','1.5')) +
        # scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels=c('0','','0.5','','1'), limits=c(0,1)) +
        labs(x=xlab, y=ylab, colour=collab) + 
        guides(colour=guide_legend(keyheight=.3, default.unit='inch'))
    p <- themefy(p)
}

Scaling function

cent <- function(v){
    v <- apply(v,2,function(x){
        x <- x - mean(unique(x))
        x <- x / max(x)
    })
}

Loading the data

allDataDir <- paste(db,'Projects/mc/data_bv3/mcBv3_maskV',sep='')
dataDirs <- dir(allDataDir)
dataDirs <- dataDirs[grep('maskV',dataDirs)]
colsOfInt <- c('participant', 'dom', 'session', 'mcBv', 'targTpeak', 'targXoff2', 'stairStart', 
               'meanRev6')
df <- data.frame()
dfRevs <- data.frame()
dfIntn <- data.frame()
for(curDir in dataDirs){
    print(curDir)
    curDf <- read.csv(paste(allDataDir,'/',curDir,'/',curDir,'.csv',sep=''))
    curDf <- curDf[,colsOfInt]
    subjStairs <- dir(paste(allDataDir,'/',curDir,'/',sep=''))
    subjStairs <- subjStairs[grep('.tsv',subjStairs)]
    for(curStairFN in subjStairs){
        curStair <- paste(allDataDir,'/',curDir,'/',curStairFN,sep='')
        curRevs <- read.table(curStair, skip=1, nrows=1)
        curIntn <- read.table(curStair, skip=4, nrows=2)
        curInfo <- readLines(curStair)
        curInfoCols <- data.frame(subj=curDf$participant[1], dom=curDf$dom[1], 
                                  sess=curDf$session[1], stairStart=curInfo[33], 
                                  mcBv=curInfo[41], targTpeak=curInfo[43],
                                  targEcc=curInfo[37], targV=curInfo[31])
        curDfRevs <- cbind(curInfoCols, curRevs[,2:9])
        nTrials <- ncol(curIntn)-1
        curDfIntn <- curInfoCols[rep(seq_len(nrow(curInfoCols)), each=nTrials),]
        curDfIntn$trial <- seq(1,nTrials)
        curDfIntn$intn <- as.numeric(curIntn[1,2:(nTrials+1)])
        curDfIntn$resp <- as.numeric(curIntn[2,2:(nTrials+1)])
        rownames(curDfIntn) <- NULL
        dfRevs <- rbind(dfRevs, curDfRevs)
        dfIntn <- rbind(dfIntn, curDfIntn)
    }
    df <- rbind(df, curDf[,colsOfInt])
}
## [1] "mc2_tgT-mcBv3_maskV-cent_p0_s1_2017-08-15_1304"
## [1] "mc2_tgT-mcBv3_maskV-cent_p11_s2_2017-08-21_1428"
## [1] "mc2_tgT-mcBv3_maskV-cent_p13_s1_2017-08-21_1755"
## [1] "mc2_tgT-mcBv3_maskV-cent_p14_s1_2017-08-22_1045"
## [1] "mc2_tgT-mcBv3_maskV-cent_p17_s1_2017-08-23_1226"
## [1] "mc2_tgT-mcBv3_maskV-cent_p18_s2_2017-08-23_1434"
## [1] "mc2_tgT-mcBv3_maskV-cent_p19_s2_2017-08-25_1037"
## [1] "mc2_tgT-mcBv3_maskV-cent_p2_s1_2017-08-16_1036"
## [1] "mc2_tgT-mcBv3_maskV-cent_p20_s1_2017-08-25_1242"
## [1] "mc2_tgT-mcBv3_maskV-cent_p23_s1_2017-08-28_1054"
## [1] "mc2_tgT-mcBv3_maskV-cent_p25_s1_2017-08-29_1349"
## [1] "mc2_tgT-mcBv3_maskV-cent_p27_s2_2017-08-30_1403"
## [1] "mc2_tgT-mcBv3_maskV-cent_p30_s2_2017-09-01_1220"
## [1] "mc2_tgT-mcBv3_maskV-cent_p5_s1_2017-08-16_1435"
## [1] "mc2_tgT-mcBv3_maskV-cent_p7_s1_2017-08-16_1637"
## [1] "mc2_tgT-mcBv3_maskV-cent_p8_s2_2017-08-22_1430"
## [1] "mc2_tgT-mcBv3_maskV-peri_p0_s2_2017-08-15_1447"
## [1] "mc2_tgT-mcBv3_maskV-peri_p11_s1_2017-08-21_1331"
## [1] "mc2_tgT-mcBv3_maskV-peri_p12_s1_2017-08-21_1544"
## [1] "mc2_tgT-mcBv3_maskV-peri_p14_s2_2017-08-22_1131"
## [1] "mc2_tgT-mcBv3_maskV-peri_p15_s1_2017-08-22_1317"
## [1] "mc2_tgT-mcBv3_maskV-peri_p17_s2_2017-08-24_1322"
## [1] "mc2_tgT-mcBv3_maskV-peri_p18_s1_2017-08-23_1348"
## [1] "mc2_tgT-mcBv3_maskV-peri_p19_s1_2017-08-24_1038"
## [1] "mc2_tgT-mcBv3_maskV-peri_p20_s2_2017-08-28_1216"
## [1] "mc2_tgT-mcBv3_maskV-peri_p23_s2_2017-08-28_1134"
## [1] "mc2_tgT-mcBv3_maskV-peri_p25_s2_2017-08-29_1431"
## [1] "mc2_tgT-mcBv3_maskV-peri_p27_s1_2017-08-30_1325"
## [1] "mc2_tgT-mcBv3_maskV-peri_p3_s1_2017-08-16_1245"
## [1] "mc2_tgT-mcBv3_maskV-peri_p30_s1_2017-09-01_1138"
## [1] "mc2_tgT-mcBv3_maskV-peri_p7_s2_2017-08-21_1630"
## [1] "mc2_tgT-mcBv3_maskV-peri_p8_s1_2017-08-16_1723"
ds <- rename(df, c(participant='subj', session='sess', meanRev6='thresh',
                   mcBv='maskV', targXoff2='targEcc'))
ds$maskV <- round(ds$maskV * 60 / 35, 1)
ds$maskV[ds$maskV<0.05] <- 0
ds$maskType <- 'fast'
ds$maskType[ds$maskV<5] <- 'slow'
ds$maskType[ds$maskV<=0.01] <- 'static'
dfIntn$maskV <- round(as.numeric(levels(dfIntn$mcBv))[dfIntn$mcBv] * 60 / 35, 1)
dfIntn$maskV[dfIntn$maskV<0.05] <- 0
ds$targEcc <- round(ds$targEcc/35,1)
ds$targType <- 'foveal'
ds$targType[ds$targEcc>1] <- 'perifoveal'
head(ds)
##   subj dom sess maskV targTpeak targEcc stairStart    thresh maskType targType
## 1    0   1    1   8.2       1.5     0.8         -2 -1.681667     fast   foveal
## 2    0   1    1  16.5       1.5     0.8          0 -1.678333     fast   foveal
## 3    0   1    1   1.0       1.5     0.8         -2 -1.448333     slow   foveal
## 4    0   1    1   0.0       1.0     0.8         -2 -1.921667   static   foveal
## 5    0   1    1   1.0       1.0     0.8          0 -1.378333     slow   foveal
## 6    0   1    1   2.1       1.0     0.8          0 -1.378333     slow   foveal

Staircase averages & convenience variables

By velocity bandwidth

thresh <- ddply(ds, .(subj,dom,sess,maskV,targTpeak,targEcc,targType), summarise,
                threshMean = mean(thresh))
thresh$maskV <- factor(thresh$maskV)

By mask type

threshType <- ddply(ds, .(subj,dom,sess,maskType,targTpeak,targEcc,targType), summarise,
                    threshMean = mean(thresh))

Centered data set

dsc <- ds
centCols <- c('dom','targTpeak','stairStart','targEcc')
dsc[,centCols] <- cent(dsc[,centCols])
dsc$maskV <- (dsc$maskV / max(dsc$maskV)) #*2 - 1
dsc$subj <- as.factor(dsc$subj)
dsc$slowC <- -1 # scaled and centered variable [-1,1]
dsc$slowC[dsc$maskType=='slow'] <- 1 #0
dsc$slowC[dsc$maskType=='fast'] <- 999
head(dsc)
##   subj dom sess      maskV targTpeak targEcc stairStart    thresh maskType targType slowC
## 1    0   1    1 0.49696970         1      -1         -1 -1.681667     fast   foveal   999
## 2    0   1    1 1.00000000         1      -1          1 -1.678333     fast   foveal   999
## 3    0   1    1 0.06060606         1      -1         -1 -1.448333     slow   foveal     1
## 4    0   1    1 0.00000000         0      -1         -1 -1.921667   static   foveal    -1
## 5    0   1    1 0.06060606         0      -1          1 -1.378333     slow   foveal     1
## 6    0   1    1 0.12727273         0      -1          1 -1.378333     slow   foveal     1

Stair plots

(p <- ggplot(dfIntn[dfIntn$targTpeak==0.5,], 
             aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ maskV, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

(p <- ggplot(dfIntn[as.numeric(as.character(dfIntn$targTpeak))==1.0,], 
             aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ maskV, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

(p <- ggplot(dfIntn[dfIntn$targTpeak==1.5,], 
             aes(x=trial, y=intn, colour=targEcc, linetype=stairStart)) +
    facet_wrap(subj ~ maskV, ncol=6) + geom_point() + geom_line() + ylim(-2,0))

Plotting across mask bandwidths

Foveal targets

ss <- sumFn(thresh[thresh$targType=='foveal',])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pCent <- themefy(plotAve(ssAve)) + ylim(yLim)
pCentNorm <- themefy(plotAve(ssAve, ylab=yLab2, yStr='norm', seStr='seNorm')) + 
    ylim(yLimNorm)

Individual plots

(plotIndiv(ssIndiv))

Perifoveal targets

ss <- sumFn(thresh[thresh$targType=='perifoveal',])
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pPeri <- themefy(plotAve(ssAve)) + ylim(yLim)
pPeriNorm <- themefy(plotAve(ssAve, ylab=yLab2, yStr='norm', seStr='seNorm')) +
    ylim(yLimNorm)

Individual plots

(plotIndiv(ssIndiv))

Group plots

w <- .56
grid.arrange(pCent + theme(legend.position='none') + ggtitle('a: foveal target'),
             pPeri + ggtitle('b: perifoveal target'),
             ncol=2, widths=c((1-w)/2,w/2))
## Warning: Removed 3 rows containing missing values (geom_linerange).

Plotting across mask types

Foveal targets

ss <- sumFn(threshType[threshType$targType=='foveal',], grpStr='maskType')
ss$maskType <- ordered(ss$maskType, levels=c('static','slow','fast'))
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pCentType <- themefy(plotAve(ssAve, grpStr='maskType', collab=colLabType)) + ylim(yLim)
pCentNormType <- themefy(plotAve(ssAve, ylab=yLab2, yStr='norm', seStr='seNorm',
                                 grpStr='maskType', collab=colLabType)) + ylim(yLimNorm)

Perifoveal targets

ss <- sumFn(threshType[threshType$targType=='perifoveal',], grpStr='maskType')
ss$maskType <- ordered(ss$maskType, levels=c('static','slow','fast'))
ssAve <- ss[ss$subj=='average',]
ssIndiv <- ss[ss$subj!='average',]
pPeriType <- themefy(plotAve(ssAve, grpStr='maskType', collab=colLabType)) + ylim(yLim)
pPeriNormType <- themefy(plotAve(ssAve, ylab=yLab2, yStr='norm', seStr='seNorm', 
                             grpStr='maskType', collab=colLabType)) + ylim(yLimNorm)

Group plots

grid.arrange(pCentType + theme(legend.position='none') + ggtitle('a: foveal target'),
             pPeriType + ggtitle('b: perifoveal target'),
             ncol=2, widths=c((1-w)/2,w/2))

Analyses

GLM

ssc <- dsc[dsc$maskType!='fast',]
pvalfn(lmer(thresh ~ dom + stairStart + sess + slowC * targTpeak * targEcc + 
                 (1|subj), data=ssc))
##                           estm    se      df    tVal  pVal star
## (Intercept)             -1.015 0.091  18.901 -11.198 0.000  ***
## dom                     -0.050 0.088  17.066  -0.567 0.578     
## stairStart               0.026 0.007 740.063   3.851 0.000  ***
## sess                    -0.010 0.015 742.864  -0.686 0.493     
## slowC                   -0.017 0.008 740.063  -2.190 0.029    *
## targTpeak               -0.037 0.010 740.063  -3.935 0.000  ***
## targEcc                 -0.051 0.008 743.272  -6.030 0.000  ***
## slowC:targTpeak          0.024 0.010 740.063   2.495 0.013    *
## slowC:targEcc           -0.040 0.008 740.063  -5.158 0.000  ***
## targTpeak:targEcc        0.000 0.010 740.063  -0.008 0.994     
## slowC:targTpeak:targEcc -0.002 0.010 740.063  -0.228 0.819

BF

Mask type

bfBase <- lmBF(thresh ~ subj + dom + stairStart + sess, data=ssc, whichRandom='subj')
bfTest <- lmBF(thresh ~ subj + dom + stairStart + sess + slowC, data=ssc, 
               whichRandom='subj')
as.vector(bfTest / bfBase)
## subj + dom + stairStart + sess + slowC 
##                               1.191478

maskTypeS X targEcc

bfBase <- lmBF(thresh ~ subj + dom + stairStart + sess + slowC + targTpeak + 
                   targEcc, data=ssc, whichRandom='subj')
bfTest <- lmBF(thresh ~ subj + dom + stairStart + sess + slowC + targTpeak + 
                   targEcc + targEcc*slowC, data=ssc, whichRandom='subj')
as.vector(bfTest / bfBase)
## subj + dom + stairStart + sess + slowC + targTpeak + targEcc + targEcc * slowC 
##                                                                       44881.41

Three masks (0,1,16)

ss <- ds[ds$maskV==0 | ds$maskV==1 | ds$maskV==16.5,]
ss$maskType <- factor(ss$maskType, c('static','slow','fast'))
ss$targType <- factor(ss$targType)
ss$targTpeak_c <- ss$targTpeak - .5
ss$targEcc_c <- (ss$targEcc - .8) / 2.1
pvalfn(lmer(thresh ~ dom + stairStart + sess + maskType * targTpeak_c * targEcc_c +
                 (1|subj), data=ss))
##                                      estm    se      df   tVal  pVal star
## (Intercept)                        -0.853 0.153  18.935 -5.574 0.000  ***
## dom                                -0.106 0.174  17.076 -0.612 0.548     
## stairStart                          0.023 0.007 544.073  3.277 0.001  ***
## sess                               -0.005 0.016 546.486 -0.326 0.744     
## maskTypeslow                        0.157 0.039 544.073  4.026 0.000  ***
## maskTypefast                       -0.276 0.039 544.073 -7.059 0.000  ***
## targTpeak_c                        -0.127 0.043 544.073 -2.956 0.003   **
## targEcc_c                          -0.028 0.040 544.651 -0.704 0.482     
## maskTypeslow:targTpeak_c            0.061 0.061 544.073  1.006 0.315     
## maskTypefast:targTpeak_c            0.056 0.061 544.073  0.933 0.351     
## maskTypeslow:targEcc_c             -0.226 0.055 544.073 -4.082 0.000  ***
## maskTypefast:targEcc_c             -0.056 0.055 544.073 -1.009 0.313     
## targTpeak_c:targEcc_c               0.008 0.061 544.073  0.139 0.890     
## maskTypeslow:targTpeak_c:targEcc_c  0.043 0.086 544.073  0.496 0.620     
## maskTypefast:targTpeak_c:targEcc_c  0.041 0.086 544.073  0.484 0.628

Bayesian test

ss$subj <- as.factor(ss$subj)
bfBase <- lmBF(thresh ~ subj + dom + stairStart + sess + maskType + targTpeak_c + 
                   targEcc_c, data=ss, whichRandom='subj')
bfTest <- lmBF(thresh ~ subj + dom + stairStart + sess + maskType + targTpeak_c + 
                   targEcc_c + maskType*targTpeak_c, data=ss, whichRandom='subj')
as.vector(bfTest / bfBase)
## subj + dom + stairStart + sess + maskType + targTpeak_c + targEcc_c + maskType * targTpeak_c 
##                                                                                    0.2685797

Interpretation: No significant difference of time-dependent slopes among mask speeds, possibly due to low power, since the analysis above finds a difference.